High Resolution Urban Heat Analysis - Correlating Measurements Between Household UAVs and Landsat Satellites

Author: Bryan Rickens

Email: rickensbj@g.cofc.edu

import pandas as pd

Part 1 - Introduction

This is a step-by-step walkthrough of my research to enhance the resolution of traditional space satellites through the application of civillian drones and remote sensing.

As remote-sensing becomes more actively utilized in the environmental sciences, our research continues these efforts in using civilian UAVs and drones for land surface temperature (LST) analysis. Given the increased spatial resolution that this technology provides as compared to standard satellite measurements, we sought to further study the urban heat island (UHI) effect -- specifically when it comes to heterogeneous and dynamic landscapes such as the Charleston peninsula. Furthermore, we sought to develop a method to enhance the spatial resolution of publicly available LST temperature data (such as those measured from the Landsat satellites) by building a machine learning model utilizing remote-sensed data from drones.

This model was constructed using data from three days (after multiple days had to be dropped due to cloud cover/mechanical error) -- 12/14/2021, 03/20/2022, and 03/28/2022.

The figure above demonstrates the problem in its simplest format: a single pixel from a Landsat satellite measurement encompasses a 900 m² area, while a civilian drone has a much more precise spatial resolution -- compared to the same Landsat area, drone images consist of many more pixels which ultimately leads to a more accurate measurement of land surface temperature.

Methodology

Echoing the problem statement, this research revolves around the attempts to enhance the resolution of Landsat satellites using civilian drones. The design of the machine learning model was to correlate the two -- by taking Landsat LST measurements and land cover (LC) classifications as inputs, and predicting the more precise drone LST readings as an output. Our end goal with our model was to be able to take any similar geographical area and to "enhance" the publicly available Landsat measurements to the more accurate readings it was trained on with drones.

Study Area

Harbor Walk - Charleston, SC

_{Taken from Google Maps}

Harbor Walk is a heterogeneous area within Charleston, South Carolina that presents a unique opportunity to study the implications of high spatial resolution in urban heat analysis. Within a short distance it consists of highly impervious city structures like concrete buildings and streets, parks with open spaces and vegetation, and the harbor itself with open water adjacent to the Cooper River. While still exhibiting the greater urban heat island effect, it is a dynamic area with many various land cover types contained within a short distance that may not be as accurately analyzed with traditional medium-to-large spatial resolutions such as 30, 60, or 100 squared meters.

Dataset Information

Landsat: Landsat LSTs were accessed from the USGS website EarthExplorer under the Collection 2 dataset.

Drone: Drone flights were conducted by Barbara Szwabowski using a DJI Mavic 2 Enterprise Dual drone.

Land Cover Classification: Land covers were accessed from the NLCD 2019 database.

Part 2 - Preprocessing

(Due to both the memory size as well as the amount of the data at hand, most preprocessing features will be explained via their respective functions)

The geospatial raster Python library, Rasterio, was used throughout the project to read and manipulate satellite/drone imagery.

A) Landsat

_{Landsat view of the Charleston Peninsula}

First, the Landsat LST needed to be updated with a scaling formula to output the values to Celsius.

import rasterio as rio

def read_landsat_image(landsat_image_path):
    """This function scales the unidentified integers in the base raster to Celsius
    input: path to Landsat file
    output: array of LSTs """
    
    with rio.open(landsat_image_path) as img:
        image = img.read(1)
        
    # converting temperatures from Kelvin to Celcius
    image_array = (image * 0.00341802 + 149.0) - 273.15

    #Outputs a temperature matrix where each pixel is the LST measurement for that area
    return image_array

Then, using the GPS coordinates of the Harbor Walk area (taken with Exiftool from the drone images to match the precise area), Rasterio was used to match the pixel locations in the temperature matrix.

from rasterio.crs import CRS
from rasterio import warp
from rasterio import transform

def match_pixel_location(landsat_image_path, long, lat, landsat_lst):
    """This function uses rasterio to correctly map the previous Landsat LST output
    to the respective geographical area required.
    
    input: 
        orginal Landsat image path, 
        list of longitudes,
        list of latitutdes,
        Landsat LST array 
        
    output:
        list of LST values for respective area """
    
    landsat = rio.open(landsat_image_path)
    dst_crs = landsat.crs #CRS from input Landsat TIF file
    src_crs = CRS.from_epsg(4326) #Standard Decimal Format for Lat/Long
    
    #Reprojects the Lat/Long coordinate value to the TIF's native format
    x,y = warp.transform(src_crs, dst_crs, long, lat)
    
    #Finds the row/col (closest center pixel location) for the value found
    row, col = transform.rowcol(landsat.transform, x, y)
    
    #Using the row/col location, find the LSTs for the given area
    lst = []
    for i in range(len(row)):
        lst.append(float(landsat_lst[row[i]:row[i]+1,col[i]:col[i]+1]))

    #Return list of Landsat LST values
    return lst

_{Cell Partition for each Landsat Pixel in the Harbor Walk Area}

B) Drone

An unorthodox method was required to take the temperature matrix given for each individual drone image to be further used. At the time of research, recent DJI models (such as the one used) required the use of third-party APIs to decode thermal metadata.

_{Drone Flight Path}

_{Thermal Image taken during flight}

def read_dji_image(image_path):
    """This function transforms a DJI Mavic 2 Enterprise Due drone image into a numpy array.
    input: JPG thermal image path.
    output: temperature numpy array."""

    # Uses the library thermal_base to process the DJI image.
    # (https://github.com/detecttechnologies/thermal_base)

    image = ThermalImage(image_path=image_path, camera_manufacturer="dji")

    thermal_np = image.thermal_np           # The temperature matrix as a np array
    raw_sensor_np = image.raw_sensor_np     # The raw thermal sensor excitation values # as a np array
    meta = image.meta                       # Image metadata

    return thermal_np, meta

The temperature matrices for each drone image were then saved as a .tiff file (where each pixel value was the LST). To be further used, Exiftool was used again to map the metadata from the original drone images to the new ones.

ArcGIS Pro was used to construct orthomosaics. Recalling the picture at the beginning, it was important to construct a singular image of each of the drone pictures stitched together to allow a one-to-one analysis between Landsat and drone LST measurements (where geogrpahical area matched exactly).

_{Constructed Orthomosaic}

After, a similar function used to find the Landsat pixels were also used to match the same area in the orthomosaics.

import numpy as np
from PIL import Image

def match_ortho_pixels(landsat_image_path, long, lat, ortho_image_path):
    """This function operates in two main parts:
        A)Taking the original Landsat, find the "corner" coordinates for each cell partition
        B)Now using the Ortho image, use these corners to find and match the cell areas -- and then output the mean LST for each
    input: 
        landsat image path,
        longitude list,
        latitude list,
        ortho image path
    output: 
        list of ortho LSTs for each cell"""
        
    landsat = rio.open(landsat_image_path)

    dst_crs = landsat.crs #CRS from Landsat TIF file
    src_crs = CRS.from_epsg(4326) #Standard Decimal Format for Lat/Long

    #Reprojects the Lat/Long coordinate value to the TIF's native format
    x,y = warp.transform(src_crs, dst_crs, long, lat)

    #Finds the row/col (closest center pixel location) for the value found
    row, col = transform.rowcol(landsat.transform, x, y)

    xUL, yUL = transform.xy(landsat.transform, row, col, offset = 'ul')
    xLR, yLR = transform.xy(landsat.transform, row, col, offset = 'lr')

    longUL, latUL = warp.transform(dst_crs, src_crs, xUL, yUL)
    longLR, latLR = warp.transform(dst_crs, src_crs, xLR, yLR)

    #Now, using the "corner" coordinates gained from the Landsat, find the corresponding Ortho cells
    ortho = rio.open(ortho_image_path)

    dst_crs = landsat.crs #CRS from ortho TIF file

    #Find the upper/left most pixel locations for each cell
    ortho_x1, ortho_y1 = warp.transform(src_crs, dst_crs, longUL, latUL)
    row_upper, col_left = transform.rowcol(ortho.transform,ortho_x1, ortho_y1)

    #Find the lower/right most pixel locations for each cell
    ortho_x2, ortho_y2 = warp.transform(src_crs, dst_crs, longLR, latLR)
    row_lower, col_right = transform.rowcol(ortho.transform,ortho_x2, ortho_y2)

    #Combined, we have the entire area for each cell.
    
    #Now convert the Ortho as an array (where each pixel is represented as its LST measurement)
    ortho_array = np.array(Image.open(ortho_image_path))
    ortho_lst = []

    #Using the area range found for each cell, find the mean (representing the mean LST for each cell area)
    for i in range(len(row)):
        ortho_lst.append(ortho_array[row_upper[i]:row_lower[i], col_left[i]:col_right[i]].mean())

    return ortho_lst

C) Land Cover Classification

Using the NLCD 2019 dataset, ArcGIS Pro was once again used to load the information and get a map of the Harbor Walk area (and its associated land covers).

The table below condenses the legend as it applies to the Harbor Walk area.

Below, you can see Harbor Walk area map as it applies to this classification legend.

_{NLCD Land Cover Classification for Harbor Walk}

Continuing to use the similar foundation used with rasterio to match the cells to the same area across all data, each cell was then further analyzed to measure its RGB proportion -- where each different color (an LULC classification) represented the cells' land cover.

def match_nlcd_pixels(landsat_image_path, long, lat, nlcd_image_path, df):
    """This function also operates in similar two main parts:
    A)Taking the original Landsat, find the "corner" coordinates for each cell partition
    B)Now using the NLCD image, use these corners to find and match the cell areas -- and then output the LC proportions for each
input: 
    landsat image path,
    longitude list,
    latitude list,
    NLCD image path,
    data frame to output to
output: 
    list of LC proportions for each cell"""

    landsat = rio.open(landsat_image_path)

    dst_crs = landsat.crs #CRS from Landsat TIF file
    src_crs = CRS.from_epsg(4326) #Standard Decimal Format for Lat/Long

    #Reprojects the Lat/Long coordinate value to the TIF's native format
    x,y = warp.transform(src_crs, dst_crs, long, lat)

    #Finds the row/col (closest center pixel location) for the value found
    row, col = transform.rowcol(landsat.transform, x, y)

    xUL, yUL = transform.xy(landsat.transform, row, col, offset = 'ul')
    xLR, yLR = transform.xy(landsat.transform, row, col, offset = 'lr')

    longUL, latUL = warp.transform(dst_crs, src_crs, xUL, yUL)
    longLR, latLR = warp.transform(dst_crs, src_crs, xLR, yLR)

    #Now, using the "corner" coordinates gained from the Landsat, find the corresponding NLCD cells
    nlcd = rio.open(nlcd_image_path)

    dst_crs = landsat.crs #CRS from NLCD TIF file
    
    #Find the upper/left most pixel locations for each cell
    nlcd_x1, nlcd_y1 = warp.transform(src_crs, dst_crs, longUL, latUL)
    row_upper, col_left = transform.rowcol(nlcd.transform,nlcd_x1, nlcd_y1)

    #Find the lower/right most pixel locations for each cell
    nlcd_x2, nlcd_y2 = warp.transform(src_crs, dst_crs, longLR, latLR)
    row_lower, col_right = transform.rowcol(nlcd.transform,nlcd_x2, nlcd_y2)

    #Combined, we have the entire area for each cell. 
    
    #Now, we convert the NLCD image to an array (where each pixel is represented by an RGB array)
    nlcd_array = np.array(Image.open(nlcd_image_path))
    nlcd_prop = []
    
    HighIntensity= np.array([171,0,0]) #Dark Red
    MediumIntensity = np.array([235,0,0]) #Light Red
    LowIntensity = np.array([217,146,130]) #Pink
    OpenWater = np.array([70,107,159]) #Blue
    OpenSpace = np.array([222,197,197]) #White
    

    for i in range(len(row)):
        sub = nlcd_array[row_upper[i]:row_lower[i], col_left[i]:col_right[i]]
        shape = sub.shape
        size = shape[0] * shape[1]

        HI = 0
        MI = 0
        LI = 0
        OS = 0
        OW = 0

        for j in range(int(shape[0])):
            for k in range(int(shape[1])):
                if (np.array_equal(sub[j,k], HighIntensity)):
                    HI += 1
                if (np.array_equal(sub[j,k], MediumIntensity)):  
                    MI += 1
                if (np.array_equal(sub[j,k], LowIntensity)):
                    LI += 1
                if (np.array_equal(sub[j,k], OpenSpace)):
                    OS += 1
                if (np.array_equal(sub[j,k], OpenWater)):
                    OW += 1

        # Output to dataframe, where each cell is updated with its measured land cover proportion
        df['HighIntensityPerc'].iloc[i] = round(HI/size,2)
        df['MediumIntensityPerc'].iloc[i] = round(MI/size,2)
        df['LowIntensityPerc'].iloc[i] = round(LI/size,2)
        df['OpenWaterPerc'].iloc[i] = round(OW/size,2)
        df['OpenSpacePerc'].iloc[i] = round(OS/size,2)


    
        
        

Unfortunately, the low spatial resolution problem is inherently baked into NLCD classifications as well -- since they also use Landsat sattelites as their basis. This can be visualized by comparing the same cell areas between drone images and the NLCD map:

_{Overlay of NLCD vs Actual Sattelite Topograpgy (via Google Earth 2022)}

As you can see, there are many areas (specifically around the borders where one LC becomes another) that do not reflect the actual geography of the area. There was even the case of a complete misclassification, where a large portion of the Charleston Aquarium was labeled as Open Space. For these reasons, multiple cells had to be dropped from analysis. You can see the dropped cells below.

_{Dropped Cells}

Part 3 - Analysis

Below you can see the overall data layout:

_{Data Overview (HI = High, MI = Medium, W = Water)}

_{Number of Cells per Set}

_{Example Data Frame Layout}

The middle date (03/20/2022) was chosen as the testing data set to avoid any errors of extrapolation. Seeing as each date had increasing temperatures, we wanted to ensure the model was trained on the highest possible LST range -- and that the testing set would contain temperatures within this range.The graph below illustrates this rationale:

_{LST Ranges for Each Date}

It is important at this point to reaffirm the limitations of this small data set: ideally, a model would be best trained on an entire year's worth of data with temperature ranges from each season. At best, this research seeks to act as a proof of concept in a limited capacity.

A) LST Correlation

Below shows a scatter plot of Landsat LST by Drone LST. While the land covers for Water and Medium Intensity exhibit a close linear relationship, High Intensity land covers have a much wider range of error when comparing measurements.

_{Landsat LST by Drone LST Correlation}

Noting this, a deeper inspection reveals some concerning differences between both Landsat and Drone LST measurements.

_{LST Distribution by Land Cover for both Landsat and Drone}

While Landsat satellites have an equally varied distribution across land cover types, Drones measurements exhibit a unique difference. While cells with Water and Medium classifications are more precisely measured with Drones, High Intensity land covers have an extreme degree of variance for LST. We believe this reaffirms our criticism of the low-spatial resolution problem -- and how a dynamic landscape such as Harbor Walk suffers from it. A cell with a 900 meter squared area cannot so easily be homogenized by land cover. We believe this extreme variance for High Intensity classifications is occuring due to the varied presence of other topography along with the impervious structures (a heavy presence of vegetation, trees, or water despite still being classified as "High Intensity" by standard 30 meter spatial resolution techniques).

This variation between LULC is especially apparent when analyzing the mean temperatures of each cell area across all days. Of particular note are the High Intensity labeled cells 8 and 18, which show significant difference with a number of other cells as seen below:

_{Error Plots of Drone LST by Each Cell}

After closer inspection, these cells are actually the only areas that solely cover a rooftop area (which tracks why they would have such a higher LST than the rest due to the highly impervious surface area).

_{Cells 8 and 18 (Roofs)}

B) Prediction Model Construction

Despite the concerns that were discovered in the analysis period, we continued forward with an attempt to use machine learning to try and find a pattern between the two forms of measurement, and ultimately attempt to "enhance" the resolution of a Landsat image by building a model used to predict the more accurate Drone LST measurement.

While various regression models were tested (such as support vector machines and tree-based/interaction methods), generalized linear models (specifically OLS) performed the best most likely due to the relatively linear relationship seen when comparing Landsat to Drone. The model was constructed and implemented using the Scikit-Learn library in Python. To see the coding process used in constructing this model, continue in this subsection. Otherwise, skip ahead to subsection C) to see the Results.

import pandas as pd

parent = pd.read_csv('ProjectData.csv')

# Seperate training/testing by date
train = parent.loc[(parent['Date'] == '14-Dec') | (parent['Date'] == '28-Mar')]
test = parent.loc[parent['Date'] == '20-Mar']

# Custom function to One Hot Encode (OHE) any categorical column in a data frame

def encode_and_bind(original_dataframe, feature_to_encode):
    """ This function one-hot encodes the desired feature column in a data frame
input: 
    original data frame,
    feature column,
output: 
    one hot encoded feature columns """
    
    # Get dummies for feature column using panda's integrated method
    dummies = pd.get_dummies(original_dataframe[[feature_to_encode]])
    # Concatenate the new OHE columns to the dataframe
    res = pd.concat([original_dataframe, dummies], axis=1)
    # Drop the original feature
    res = res.drop([feature_to_encode], axis=1)
    # Return the updated data frame
    return(res) 

(OHE was utilized to avoid any unintentional linear relationships commonly seen through Label Encoding)

# Apply OHE function to both sets
train = encode_and_bind(train, 'LULC')
test = encode_and_bind(test,'LULC')

# To avoid multicollinearity/Dummy Variable Trap, drop one category
train.drop(columns='LULC_Medium', inplace = True)
test.drop(columns='LULC_Medium', inplace = True)

# Now, assign between X and y train/test variables
X_train = train[['LandsatTemps','LULC_High','LULC_Water']]
y_train = train['DroneTemps']

X_test = test[['LandsatTemps','LULC_High','LULC_Water']]
y_test = test['DroneTemps']

# OLS model through sklearn

from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error

model = LinearRegression()

# Fit on training
model.fit(X_train,y_train)

# Predict on test
y_pred = model.predict(X_test)

# Find RMSE
rmse = mean_squared_error(y_test,y_pred,squared = False)

For a more robust model analysis, a results data frame was also created using aspects of the original master data frame. This was further utilized in the next Results section.

# Construct initial data frame
results = pd.DataFrame( {
    'Actual Drone LST' : y_test,
    'Predicted Drone LST' : y_pred,
    'Difference' : y_test - y_pred
})

# Concat useful information from parent data frame
results['CellNum'] = parent['CellNum']
results['Date'] = parent['Date']
results['LULC'] = parent['LULC']
results['Landsat LST'] = parent['LandsatTemps']

C) Prediction Model Results

The linear model resulted in an explained variance of 49% with an RMSE of 1.01 degrees Celsius.

Using Shapley values (with the SHAP library in Python) to assess feature importance, it is clear that Landsat Temp had the dominating influence on our linear model. In fact, removing LULC as a variable only resulted in a marginal difference (with an R^2 of 0.42 and an RMSE of 1.07).

_{Shapley Values of Variables}

Following what was seen previously, residuals of cells of Open Water and Medium Intensity were much smaller than those of High Intensity.

_{RMSE by Land Cover Classificaton}

With the validity of the model in question, the following figures will demonstrate key errors undermining model performance. Below demonstrates the Actual vs Predicted values, and we see that High Intensity predictions follow no linear trend and exhibit high margins of error.

_{Actual vs Predicted Values}

Similarly, the next figure highlights the heteroscedasticity of High Intensity predictions. While Water and Medium predictions are more precise, predictions of High land covers vary wildly from the actual values.

_{Residuals of Predicted Values}

Visualizing the residuals by the actual values as well further illuminates the problems with the linear model when it comes to High Intensity land cover predictions. Interestingly, the model seems to severely undervalue all High Intensity predictions up to a certain threshold (29 degrees Celsius). Immediately past this threshold, the model inversely begins to overvalue the predictions.

_{Residuals of Actual Values}

This could also be further evidence of the low spatial resolution problem bearing its head again. As you can see, Water and Medium LST have distinctly different ranges that do not overlap. High Intensity, however, exhibits LST values across the entire distribution: with some cells close in temperature to Water cells, some similar to Medium cells, and others higher than the rest of the cells. It seems that for the High Intensity cells that fall in the range of other land cover types, the model under predicts; for cells within their own high temperature range, the model then over predicts.

To further investigate what we assumed to be a consequence of low spatial Landsat resolution, we identified the cells with the highest residuals and visually examined the geography of the area in comparison to their official NLCD land cover classification. Below are three cells from each land cover type that exhibited the highest residuals in their group (excluding the two roof cells for High that were previously discusssed).

_{Cell 23 (High Intensity)}

_{Cell 24 (Medium Intensity)}

_{Cell 15 (Water)}

Each group example exhibits the same issue: the bordering presence of competing topography of another land cover. Because of its low spatial resolution and its wide 900 square meter area coverage, traditional Landsat land cover classification methods used on heterogeneous landscapes such as Harbor Walk are seen to be unfit. Because of the landscape's highly dynamic nature that includes various land cover types over a small distance, traditional methods homogenize the area which fails to capture their true nature.

• Cell 23 (High Intensity) includes heavy vegetation such as trees and the outskirts of a park
• Cell 24 (Medium Intensity) includes a street in the upper corner amongst the vegetation and park area
• Cell 15 (Water) includes a significant portion of a building at the edge of the open water

While potentially still following the metrics outlined in the NLCD legend listed earlier (in relation to the percentage of impervious structures included), it has still been shown that especially dynamic areas can result in an increased variance of topography even within the same land cover class. This variance is not as easily measured with just one pixel's worth of area coverage as seen in traditional low spatial resolution methods such as Landsat.

Part 4 - Conclusion

We attribute a fundamental design flaw to our unsuccessful attempts: essentially using the problem in an effort to solve the problem. As discussed previously, using NLCD classifications of land covers using the very same medium spatial resolutions used in Landsat measurements proved to only further highlight the disparity between satellite and drone resolutions. While land covers classified with open water and significant vegetation can be accurately enhanced, those classified with high degrees of impervious surfaces cannot. While ineffective for our model, we believe these results underline the repercussions of contemporary LULC classifications based on medium spatial resolutions.

As seen in our results, highly dynamic and heterogeneous areas such as those in the Charleston peninsula exhibit varying correlation with publicly available LST measurements taken from a standard medium resolution. While areas with a uniform distribution of surface area (be it water or vegetation) exhibited a stronger linear relationship between drone and Landsat measurements, dynamic areas consisting of a range of land cover types within the same distance exposed increased variation between medium and high resolution measurements. Areas containing a large number of impervious surfaces have been shown to contain the most error when compared to higher resolution measurements captured with drones.

To better analyze urban heat islands and their implications on the environment, government and other publicly available evaluations of land surface temperature need to adapt to higher resolutions with a smaller spatial areas per-pixel. Given our results indicating increased disparity specifically with areas containing mostly impervious structures, urban environments are uniquely error prone using traditional land surface temperature estimates. While we tried unsuccessfully to enhance this resolution in the meantime with predictive modeling techniques, we hope to have shed further light on the issue and demonstrated the implications of the problem at large.